Non-Commutativity of the Zero Chemical Potential Limit 
and the Thermodynamic Limit in Finite Density Systems 
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Monte Carlo simulations of finite density systems are often plagued by the complex action problem. 
We point out that there exists certain non-commutativity in the zero chemical potential limit and 
the thermodynamic limit when one tries to study such systems by reweighting techniques. This is 
demonstrated by explicit calculations in a Random Matrix Theory, which is thought to be a simple 
qualitative model for finite density QCD. The factorization method allows us to understand how the 
non-commutativity, which appears at the intermediate steps, cancels in the end results for physical 
observables. 



PACS numbers: 05.10.Ln, ll.25.-w, 11.25.Sq 

I. INTRODUCTION 

QCD at finite baryon density and/or finite tempera- 
ture has attracted much attention due to its relevance 
to the physics of the early universe, heavy ion colli- 
sions and neutron stars. It is of particular importance 
to explore the phase diagram in the /i (chemical poten- 
tial), T (temperature) plane, where interesting phases 
such as a superconducting phase have been conjectured 
to appear. Monte Carlo simulations, which enable first- 
principle studies at /i = 0, are hindered by the fact that 
the fermion determinant becomes complex for fi ^ 0. 
The standard reweighting method uses the absolute value 
for generating configurations, and takes account of the 
phase in measuring observables. Due to cancellations 
caused by the oscillating phase, however, the required 
number of configurations grows exponentially with the 
system size. Furthermore the method suffers from the so- 
called overlap problem, which is caused by the mismatch 
of the region in the configuration space which contributes 
to the ensemble average and the region which one mostly 
samples in the phase-quenched model. 

One of the recent developments in this direction is that 
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the reweighting techniques have proven to be of use in ex- 
ploring the phase diagram at small fi and large T, where 
the fiuctuation of the phase is still under control. Sub- 
stantial improvements are achieved by various tricks such 
as the multi-parameter reweighting [l|, I3] , Taylor expan- 
sion approach Q and the imaginary n approach 0, |^ . 
In particular the first approach was able to locate the 
critical end point in the n-T plane (See also Ref. 0.) 

In Ref. two of the authors (K.N. A. and J.N.) pro- 
posed a new method, the factorization method, for sim- 
ulating systems with complex actions. The method uti- 
lizes the fact that the distribution of observables factor- 
izes into the corresponding distribution for the phase- 
quenched model and the weight factor representing the 
effect of the phase. Each factor can be obtained by con- 
strained Monte Carlo simulations, which eliminate the 
overlap problem completely. The knowledge of the weight 
factor allows us to understand the effect of the phase in- 
tuitively. The method is quite general, and in particular 
it is expected to be useful in going beyond the small fi 
regime in finite density QCD. The method proposed for 
simulating 0-vacuum like systems j8| may be viewed as 
a particular case of the factorization method. In Ref. 
it was pointed out that the factorization method belongs 
to the class of methods known as 'the density of states 
method'. 

In Ref. we have tested the method in a Ran- 
dom Matrix Theory (RMT), which is thought of as a 
schematic model for QCD at finite baryon density [lT| . 
RMT was originally introduced to describe the spectrum 
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of the Dirac operator and has been extensively studied 
in the hterature The model we have studied ex- 

hibits a first order phase transition at some value of the 
"chemical potential" and can be solved analytically even 
for finite size matrices One hopes to capture the 

essential properties of realQCD, while at the same time 
having a testing ground 14(1 for methods to be applied 
to real QCD. Indeed the factorization method is quite 
successful in RMT where exact results are repro- 
duced with great accuracy, and one understands clearly 
how the phase of the determinant induces the first order 
phase transition. 

In this paper we investigate this model further in or- 
der to address the question of the non-commutativity in 
the /i — > limit and the thermodynamic limit. The 
factorization method allows us to understand how the 
non-commutativity appears at the intermediate steps of 
reweighting techniques, and how it cancels in the end re- 
sults for physical observables. Preliminary results have 
already been presented in conference proceedings, Refs. 

Cilia. 



II. RMT FOR FINITE DENSITY QCD AND 
THE MONTE CARLO METHODS 

We consider the RMT defined by the partition function 

Z = J dH^e-^*'-^^''^) detD , (II. 1) 

where is a x complex matrix, and _D is a 2N x 2N 
matrix given by 
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The above model can be thought of as a schematic model 
of finite density QCD with one flavor, where the param- 
eters /i and m correspond to the chemical potential and 
the quark mass respectively. The size of the matrix W 
corresponds to the number of low lying modes of the 
Dirac operator and if the density of these modes is taken 
to be unity, N can be interpreted as the volume of space- 
time. As physical observables, one may consider the "chi- 
ral condensate" and the "quark number density" defined 
by 
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tr {^iD-^) , 
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In what follows we consider the massless case (m 
and focus on the "quark number density" . 

The model was first solved in the large- A'' limit [Tll |. 
but an analytic solution has been obtained in C3l even 
for finite A''. The partition function can be expressed as 



Z{ii) 



(-1) 



7V+1 



-7(Ar + l,K) 



(II.5) 



where k = —Np? and 7(71, x) is the incomplete 7- 
function defined by 



l{n,x) 



e"* dt 



(11.6) 



From this one obtains the vacuum expectation value 
(VEV) of the quark number density as 
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iv) = iilnZ(u) 
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In Fig. n we plot {v) as a function of the chemical po- 
tential ixiov N = 8, 16, 32, 64, 128. The large- A^ hmit of 
this formula is easily found by applying the saddle-point 
method to the incomplete 7-function. We obtain 



lim 

N^oc 



—fi for /i < /ic 
l//i for 11 > fic 



(11.9) 



where /ic is the solution to the equation 1 + +ln(/i^) = 
0, and its numerical value is given by /ic = 0.527 • • •. We 
find that the quark number density (v) has a disconti- 
nuity at /i = /ic. Thus the schematic model reproduces 
qualitatively the first order phase transition expected to 
occur in 'real' QCD at nonzero baryon density. 




FIG. 1: The exact result 111.811 for the 'quark number density' 
(i/) is plotted as a function of the 'chemical potential' /i for 
iV = 8, 16, 32, 64, 128. In the A'' ^ cx) limit, the function 
develops a discontinuity ai jj, — /ic = 0.527 . ■ ■. 



The model (III. Ill cannot be simulated directly because 
the fcrmion determinant has a non-zero phase F defined 
by 



detZ? = e'^ldetDl 



(11.10) 



when fi =/= 0. In order to reveal the importance of the 
phase let us consider the phase-quenched model 



Zn = dW 



detD I 
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= J dWe-^" , 
5*0 = N tr (W^W) -ln\detD\ 



(11.11) 
(11.12) 



and the corresponding VEVs, which we denote as (. . .)o. 
The large- hmit of {i/)o can be obtained analytically as 



lim (z/)o = 



/i for fi < 1 
l/fi for /i > 1 



(11.13) 



Comparing this with the result ljll.9(l for the full model, 
we see that the effect of F is dramatic for fj, < 1. 
The reweighting method uses the identity 



(11.14) 



to calculate {ly) for the full model by Monte Carlo simula- 
tion of the phase-quenched model. Due to the symmetry 
of the phase-quenched model under W i-^ —W, where 
the fermion determinant det D as well as the observable 
becomes complex conjugate, we obtain 



(z^RCOsr)o 
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(cosr)o 

(z/isinr)^ 

(cosr)o 



(11.15) 
(11.16) 
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where ^'r and vi denote the real part and the imaginary 
part of respectively. The same symmetry implies that 



(11.18) 



For the unquenched model, on the other hand, one ob- 
tains 

(i'r)=0 ; {iA)^ifi (11.19) 

in the large N limit for fx < fi^, which is in sharp contrast 
to the phase-quenched result IjlLlSI) . 

The factorization method in the present case amounts 
to calculating the VEVs on the r.h.s. of ljII.16|) and ljII.17|) 
by the formulae 



(^'Rcosr)^ 

(cosr)o 

(i/isinr)^ 



dxxp'^\x)wR{x) , (11.20) 
dxp'^\x)wK{x) , (11.21) 
dx X p["\x) wi{x) , (11.22) 



where the functions p^f^ (x) {i — R, I) represent the dis- 
tribution of i^i [i = R, I) in the phase-quenched model 

pf\x)^^ {5{x-v,))o . (11.23) 



The weight factors Wi{x) in Eqs. (|II.20p - HII.22|l are de- 
fined by 



wi{x) 



def 



dcf 



(cOSr)R,:r , 



(sinT) 



l,X 1 



(11.24) 
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where the VEV (. . .)i^x is taken with respect to yet an- 
other partition function 



Z,{x) = j dPFe~^° 5{x - Vi) 



(11.26) 



The problem then reduces to the calculation of the four 
functions pf^\x) and Wi(x) (i = R, I). This may be done 
in various ways. In Ref. as well as in the present 
work we have simulated 



(11.27) 



where the (5-function in (|II.26|) is replaced by a sharply 
peaked potential V{x), which is taken to be Gaussian 



0' 



(11.28) 



with 7 and ^ being real parameters. By choosing 7 large 
enough, the results become insensitive to its value (we 

used 7 = 1000.0). The functions pf^\x) can be obtained 
from the same simulation. We refer the reader to Ref. 
[lol | for more details. 

The complex action problem occurs in the reweighting 
method because the trigonometric functions in ljII.16|) 
and 



(|II.17I) flip their signs violently when one moves 
around the configuration space. The same problem oc- 
curs in the factorization method when one calculates the 
weight factors ljll.24|l . Ijll.25p . However, by simulating 

(III.27II. one forces the simula- 



the constrained system 
tion to sample the important region of the configuration 
space, which is rarely visited by the simulation of the 
phase-quenched model IjlLlip . Thus the factorization 
method removes the overlap problem completely. Once 
the weight factors are obtained roughly, one can make 
the sampling more efficient by the use of multi-canonical 
simulations. This is not yet done, however. The knowl- 
edge of the weight factors is also useful in understanding 
the effect of the phase intuitively 0, 0| , and it plays a 
crucial role in the present work. 



III. NON-COMMUTATIVITY OF THE ^ ^ 
AND THE THERMODYNAMIC LIMIT 

In this Section we discuss the non-commutativity of 
the two limits /i — > and ^ 00 in the RMT. Such 
non-commutativity can be readily seen from the parti- 
tion function (|II.5p . Below the critical point, the large N 
behavior is given by Z{ij) ^ const. e~^'^ . Omitting the 
/i-independent factor, the partition function approaches 
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unity if one takes the fj, limit first, whereas it van- 
ishes if one takes the large N limit first. More generally, 
one obtains e"'^, if one takes the two limits simultane- 
ously with fixed N^'^ = C . This non-commutativity is 
caused by the phase of the determinant. The phase van- 
ishes at /i = for finite iV, and one obtains a nonzero 
result for the partition function in the large N limit with 
appropriate normalization. If one takes the large N limit 
first for small but finite ^, however, the oscillation of the 
phase becomes violent, and as a result one obtains Z = Q 
with the same normalization. 

We note, however, that the free energy defined by 

lim <! 4 In Z{fi, N) - const. I (III.29) 



N' 



after subtracting the /i-independent constant is give n b; 



n py 



fifj-) — fi at fi < iif., which is continuous at /i = 
Furthermore we find from ljll.8(l or Fig. that the quark 
number density does not have such non-commutativity, 
either. Thus the non-commutativity does not seem to 
appear in physical quantities, but it appears at the inter- 
mediate steps of the reweighting method as we will see 
in what follows. 

Let us first consider the weight factor wb,{x), which has 
a non-commutativity similar to the partition function. 
Since the phase F disappears identically for fi = 0, one 
obtains wr(x) = 1 at = for any N. On the other 
hand, one obtains w-[{{x) = in the large N limit for any 
/J, ^ 0, since the phase oscillates violently. In Fig. [21 we 
plot wk{x) for ^j, = 0.1 and /x = 0.2 at iV = 8,16,32. 
The weight factor w^,{x) crosses zero, and the crossing 
point moves to infinity as —>■ 0. Thus the convergence 
of w^{x) to WYi{x) = 1 in the fi ^ limit is not uniform. 




where the coefficient a-ii{N) grows linearly with A''. 
Therefore the asymptotic behavior of the weight factor 
w^{x) depends on how the two limits fj, 0, N ^ oo 
are taken. 
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FIG. 3: The slope a-R{f^, N) of the weight factor wn(x) in 
the linear regime is plotted against /i for N — 8, 32, 64. For 
sufficiently small /i it fits well to aR{fi,N) ~ — QfR(A'') /i. 



Let us next turn to p'^\x). In Fig. 0] we plot it for 
various iV at /i = 0.2. At small N the distribution is 
peaked near the origin and the dependence on N is small. 
As we go to larger N the peak moves to a; ^ 0.2 and starts 
to grow. The VEV (^'r)o, which represents the position 
of the peak, is therefore close to zero at small N, and it 
approaches (i^r)o = M i'^ the large N limit. 
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FIG. 2: The weight factor wn{x) is plotted for /i = 0.1 and 
0.2 at AT = 8, 16,32. 
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FIG. 4: The distribution p^^\x) of j^r in the phase-quenched 
model is plotted for /i — 0.2 at various TV. 



As /i approaches zero for fixed iV, a linear regime devel- 
ops in the region where wr(x) crosses zero. We extract 
the slope (TR(/i,iV) in this regime and plot it against fi 
in Fig. 131 At small fi we find that the slope can be fitted 
nicely by 



(III.30) 



In Fig. Owe plot (i^r)o — M against 1/N for various fi. 
We see that the large N asymptotic behavior is given by 



1 

N 



(III.31) 



where the coefficient C(/i) grows as ~ for /i 
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From this one finds that the large N scahng sets in at 

0.2 



TV » iVt, ~ — . 



(III.32) 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 
1/N 

FIG. 5: We plot {u-r}o — fJ- against 1/N for various ^. 



The product p^\x) WR(a;) gives the unnormahzed dis- 
tribution of i^R in the full model, which we plot in Fig.El 
The distribution itself, even after appropriate normaliza- 
tion, has the non-commutativity which is inherited from 
p^^\x) and wr{x). However, the VEV (i^r), which is 
the first moment of the distribution, is always close to 
zero as one can see from Table The reason depends on 
whether p, ^ or p In the former case the 

distribution is peaked around the origin, which makes 
the first moment close to zero. In the latter case the 
positive and negative regions of the distribution cancel 
each other in the calculation of the first moment. Thus 
the non-commutativity cancels in the end result for the 
VEV (i'r). The situation for the imaginary part (vi) is 
discussed in the Appendix. 
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W) 


8 


0.0056(6) 


-0.1970(5) 


-0.1915(7) 


16 


0.0060(4) 


-0.1905(13) 


-0.1845(13) 


24 


0.0076(9) 


-0.1972(14) 


-0.1896(17) 


32 


0.0021(8) 


-0.1947(19) 


-0.1927(25) 


48 


0.0086(37) 


-0.2086(54) 


-0.2000(88) 



TABLE I: The VEVs (^r), i (ui) and {v) obtained by the 
factorization method for /i = 0.2 at various N. Statistical 
errors computed by the jackknife method are shown. 
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FIG. 6: The product p^\x)w-R,{x) , which gives the unnor- 
mahzed distribution for un, in the full model, is plotted for 
p = 0.2 at various N. 



be understood as a property of the phase-quenched par- 
tition function 



II. lip . This model undergoes a phase 
transition to a phase of condensed Goldstone bosons with 
nonzero baryon number at p = mT^jl [ill Il7|. For zero 
quark mass this phase transition takes place at = 0. If 
we take the thermodynamic limit before the /i ^ limit, 
the observables are calculated in a Bose-condensed phase. 
If the limits are taken in the opposite order, the ground 
state is the regular chirally broken vacuum state without 
Bose condensate. For this reason we obtain different val- 
ues for observables depending on the order of the limits. 
For example, if the chiral condensate (E)o is calculated 
before the /i — > limit, we find (S)q = 0, whereas if we 
take the limit first, we find (E)o / 0. Another 

example is -§^{i-')q, which vanishes if the — > limit is 
taken before the thermodynamic limit, while it is equal 
to unity if the limits are taken in the opposite order. 

These two domains are separated by the weak non- 
hermiticity limit 18], where the thermodynamic limit 
is taken at fixed p^N . This is the domain where the 
real part of the eigenvalues is of the order of the aver- 
age spacing of the eigenvalues (the eigenvalues are purely 
imaginary for /x = 0), and quantities such as the average 
spectral density of the Dirac operator can be obtained 
analytically in this limit [T^ . 

Mathematically, the boundary p^N gives the "aver- 
age" radius of convergence of the perturbative expansion 
of the fermion determinant in powers of p. According 
to Kato's criterion the perturbative series is convergent 
if the norm of the perturbative operator is less than the 
level spacing. Indeed the norm of /i74 is p"^ and the aver- 
age level spacing of random matrix Dirac operator ljll.2|) 
is ~ l/N. 

The Bose condensed phase is not present in the full 
partition function (|II.1|I . Observables depend smoothly 
on p for p < pc 0. Apparently the phase oscillations of 
the fermion determinant completely wipe out this phase 
transition. However, because reweighting methods are 
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const. 
7^ 



based on the phase-quenched partition function, we see 
traces of the non-commutativity of the /i ^ hmit and 
the thermodynamic limit at intermediate steps of the cal- 
culation. 

It is expected that such non-commutativity appears 
also when one studies real QCD at finite baryon density 
by reweighting type methods, and the transition occurs 
when the system size becomes larger than Vtr 

The system size in the recent works at small /x 

01 may be below the transition point, but the transition 
will occur if one goes to larger /i for the same system size. 

Since the non-commutativity cancels in the end results 
for physical observables in the full model, the transition 
is not a physical one, but it should rather be considered 
as a property of the reweighting type methods. We hope 
that our results will be useful when one tries to go beyond 
the small ^ regime in real QCD. It should also be men- 
tioned that the conjectured superconducting phase may 
be easier to access from the other extreme, namely from 
the large /i regime, where the fluctuation of the phase 
becomes milder again according to the results in RMT 

M. 



Acknowledgments 

We are grateful to Misha Stephanov for valuable com- 
ments. J. A. acknowledges the support by the EU net- 
work on "Discrete Random Geometry" , grant HPRN- 
CT-1999-00161 as well as by MaPhySto, Network of 
Mathematical Physics and Stochastics, funded by a grant 
from Danish National Research Foundation. K.N.A.'s 
research was partially supported by RTN grants HPRN- 
CT-2000-00122, HPRN-CT-2000-00131 and HPRN-CT- 
1999-00161 and the INTAS contract N 99 0590. The work 
of J.N. was supported in part by Grant-in- Aid for Scien- 
tific Research (No. 14740163) from the Ministry of Ed- 
ucation, Culture, Sports, Science and Technology. J.V. 
was supported in part by U.S. DOE Grant No. DE-FG- 
88ER40388. 



V. APPENDIX 

In this Appendix we briefly describe the situation with 
the imaginary part vi of the quark number density. The 
weight factor wi{x) becomes wi{x) = at /x = for any 
N, and it also becomes wi{x) = in the large A'^ limit 
for any /i 7^ 0. However the two extreme cases are not 
connected smoothly. In Fig.[7|we plot the weight factor 
wi{x) for fi = 0.1 and fi = 0.2 at = 8, 16, 32. 

Since wi{x) is an odd function due to symmetry, it 
crosses the origin. A linear regime is seen to extend from 
the origin as /i goes to zero for fixed A^. We extract the 
slope in the linear regime, and plot it as a function of fi 
in Fig. |S1 At small /i, the slope can be fitted nicely by 



1 

0.9 
0.8 

0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 

-0.1 



8 0.1 ' — - 

16 0.1 ' - 

32 0.1 ■ ■ 

8 0.2 " 
16 0.2 

32 0.2 - - 



0.2 



0.4 0.6 0.8 

X 



1 



1.2 



1.4 



FIG. 7: The weight factor wi{x) is plotted for p = 0.1 and 
0.2 a,t N = 8, 16,32. 
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FIG. 8: The slope ai{fi, N) of the weight factor wi{x) in 
the linear regime is plotted against fi for N — 8, 32, 64. For 
sufficiently small n, it fits well to cri(/^, N) ~ —ai{N) fi. 
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FIG. 9: The function p[°\x) is plotted for N = 8, 16, 24, 32, 
48, 64, 96 at ^ = 0.2. 



ai{fi,N)^-ai{N)fi , (V.33) 
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where the coefScient ai{N) grows hnearly with N. Thus 
we find that the weight factor wi{x) has similar non- 
commutativity as wb.{x). 

On the other hand, p[^\x) does not depend much on 
/i, and the peak at a; = grows smoothly with N as one 



can see from Fig. |5| The end result for i {vi) does not 
have the non-commutativity (See Table 0, but the can- 
cellation in this case occurs between the numerator and 
the denominator of ljII.17|l . which makes it less obvious 
than the situation with (i'r). 
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